

#include "sl_math.h"

#include<math.h>

#define Meisen 2147483647  //梅森素数（形式为2^(p)-1的素数）：m=2^(31)-1
#define Multiplier 16807       //乘子 a=7^(5)

static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。


/// <summary>
/// 计算sin函数。x = 0 处用泰勒展开： sin x = x - (1/3!)*x^3+ (1/5!)*x^5
/// 
/// sinf(x)在五次多项式空间上的正交投影就是对sin(x)的最佳五次多项式逼近，这个近似在整个周期内的精度都非常高。
/// 这个五次多项式的表达式是  y = 0.987862x - 0.155271x ^ 3 + 0.00564312x ^ 5
/// </summary>
/// <param name="theta">弧度制角度。3.14表示180度</param>
/// <returns>sin的值</returns>
real_t c_sin(real_t theta)
{
	// sin x = x - (1/3!)*x^3+ (1/5!)*x^5
	//sinf(x)在五次多项式空间上的正交投影就是对sin(x)的最佳五次多项式逼近，这个近似在整个周期内的精度都非常高。
	//这个五次多项式的表达式是  y = 0.987862x - 0.155271x ^ 3 + 0.00564312x ^ 5

	real_t x2, x3, x5, y = 0;
	bool flag = true;  //符号位。
	if (theta < 0)
	{
		flag = !flag;
		theta = -theta;
	}
	real_t x;
	int n = (int)(theta / PI_mul_2);
	x = theta - n * PI_mul_2;
	if (x > PI) {
		flag = !flag;
		x = x - PI;
	}

	if (x == PI || x == 0)y = 0;
	else if (x < PI_div_2) {
		x2 = x * x;
		x3 = x * x2;
		x5 = x2 * x3;
		y = 0.987862f * x - 0.155271f * x3 + 0.00564312f * x5;
	}
	else if (x == PI_div_2) y = 1;
	else if (x > PI_div_2 && x < PI) {
		x = PI - x;
		x2 = x * x;
		x3 = x * x2;
		x5 = x2 * x3;
		y = 0.987862f * x - 0.155271f * x3 + 0.00564312f * x5;
	}

	if (flag == true)return y;
	else return -y;
}

/// <summary>
/// 计算cos函数。直接调用c_cos函数计算。
/// </summary>
/// <param name="theta">弧度制角度。3.14表示180度</param>
/// <returns>cos的值</returns>
real_t c_cos(real_t theta)
{
	return c_sin(theta + PI_div_2);
}


/// <summary>
/// 计算tan的反函数。
/// </summary>
/// <param name="tan">tan 函数的值</param>
/// <returns>tan对应的角度，-PI/2 ~ PI/2</returns>
real_t c_atan(real_t tan)
{
	if (ABS(tan) <= 1) {
		real_t x = tan;
		real_t ax = ABS(x);
		return PI_div_4 * x - x * (ax - 1) * (0.2447 + 0.0663 * ax);
	}
	else
	{
		real_t x = 1 / tan;
		real_t ax = ABS(x);
		real_t t= PI_div_4 * x - x * (ax - 1) * (0.2447 + 0.0663 * ax);
		if (tan > 0) {
			return PI_div_2 - t;
		}
		else
		{
			return -PI_div_2 - t;
		}
	}
	/*
	首先，考虑圆的对称性，在0-PI/4的范围内计算准确三角函数即可。其他情况可以用对称性计算。
	在0-PI/4 的区间，|tan|<=1,令x=tan，用线性近似：arctan(x)=PI/4*x. 最大误差有0.07rad = 4度。
	ax(x−1)形式的二阶近似：arctan(x)=PI/4*x+0.273x*(1−|x|),最大近似误差 0.0038 rad=0.22度
	x(x−1)(ax−b)形式的三阶近似：arctan(x)=PI/4*x−x*(|x| −1)*(0.2447 + 0.0663*|x| ),最大近似误差 0.0015 rad = 0.086度
	*/
}


/// <summary>
/// tan = y/x.  计算这个点与x轴的角度
/// </summary>
/// <param name="y">y坐标</param>
/// <param name="x">x坐标</param>
/// <returns>tan对应的角度，-PI ~ PI</returns>
real_t c_atan2(real_t y, real_t x)
{
	if (x > 0) {
		real_t tan = y / x;
		return c_atan(tan);
	}
	if (x < 0) {
		real_t tan = -y / x;
		return PI - c_atan(tan);
	}
	//if (x == 0) 
	else{
		if (y > 0)return PI_div_2;
		if (y < 0)return -PI_div_2;
		else { return 0; }
		// if (y == 0)return 0;
	}
}


/// <summary>
/// 返回均匀分布的随机数。最小标准随机数生成器
/// </summary>
/// <param name="min">随机数的最小值</param>
/// <param name="max">随机数的最大值</param>
/// <returns>均匀分布的随机数</returns>
real_t random_ave(real_t min, real_t max)
{
	//声明一个静态局部变量，这个变量实际地址不在函数中，程序编译时旧已经赋值了，调用这个程序的时候不会重复赋值。
	//static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。

	x_i = (Multiplier * x_i) % Meisen;
	return min + (max - min) * ((real_t)x_i / Meisen);  //
}


/// <summary>
/// 返回均匀分布的随机数。最小标准随机数生成器。在【-1，1】之间分布。
/// </summary>
/// <returns>均匀分布的随机数</returns>
real_t random11(void)
{
	//声明一个静态局部变量，这个变量实际地址不在函数中，程序编译时旧已经赋值了，调用这个程序的时候不会重复赋值。
	//static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。

	x_i = (Multiplier * x_i) % Meisen;
	return 2 * ((real_t)x_i / Meisen) - 1;  //
}

/// <summary>
/// 返回均匀分布的随机数。最小标准随机数生成器。在【0，1】之间分布。
/// </summary>
/// <returns>均匀分布的随机数</returns>
real_t random01(void)
{
	//声明一个静态局部变量，这个变量实际地址不在函数中，程序编译时旧已经赋值了，调用这个程序的时候不会重复赋值。
	//static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。

	x_i = (Multiplier * x_i) % Meisen;
	return  ((real_t)x_i / Meisen) ;  //
}

/// <summary>
/// 返回无符号的随机整数
/// </summary>
/// <returns></returns>
unsigned int randomint(void) {
	x_i = (Multiplier * x_i) % Meisen;
	return x_i;
}


/// <summary>
/// 高斯随机变量生成函数。使用Box_Muller方法，避免三角函数运算，速度稍微快一些
/// </summary>
/// <param name="ave">高斯分布均值</param>
/// <param name="sigma">高斯分布的标准差</param>
/// <returns>高斯随机数</returns>
real_t random_gauss_as(real_t ave, real_t sigma)
{
	//声明静态局部变量，这个变量实际地址不在函数中，程序编译时旧已经赋值了，调用这个程序的时候不会重复赋值。
	static unsigned int call_times = 0;//随机数调用次数。
	static real_t random1 = 0;
	static real_t random2 = 0; //0-1分布的随机数
	static real_t gauss_1 = 0;
	static real_t gauss_2 = 0; //均值0，方差1的高斯分布。 N(u,sigma^2)= N(0,1)
	//static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。

	call_times++;
	if (call_times % 2 == 1) {
		real_t r = 10;
		while (r > 1) // Box_Muller方法是一种拒绝方法。r如果大于1就淘汰
		{
			x_i = (Multiplier * x_i) % Meisen;
			random1 = ((real_t)x_i / Meisen);  //现在d是0-1的分布。最小标准随机数生成器
			x_i = (Multiplier * x_i) % Meisen;
			random2 = ((real_t)x_i / Meisen);  //现在d是0-1的分布

			random1 = random1 * 2 - 1;   //将0 - 1分布 变成-1 - 1的分布
			random2 = random2 * 2 - 1;


			r = random1 * random1 + random2 * random2; //刷新r
		}
		r = sqrtf(-2 * logf(r) / r);
		gauss_1 = r * random1;
		gauss_2 = r * random2;    //一次可以计算两个高斯随机变量

		return gauss_1 * sigma + ave;
	}
	else
	{
		return gauss_2 * sigma + ave;
	}
}


/// <summary>
/// 高斯随机变量生成函数。这里需要用到三角函数，速度稍慢一些。
/// </summary>
/// <param name="ave">高斯分布均值</param>
/// <param name="sigma">高斯分布的标准差</param>
/// <returns>高斯随机数</returns>
real_t random_gauss_ss(real_t ave, real_t sigma)
{
	//声明静态局部变量，这个变量实际地址不在函数中，程序编译时就已经赋值了，调用这个程序的时候不会重复赋值。
	static unsigned int call_times = 0;//随机数调用次数。
	static real_t random1 = 0;
	static real_t random2 = 0; //0-1分布的随机数
	static real_t gauss_1 = 0;
	static real_t gauss_2 = 0; //均值0，方差1的高斯分布。 N(u,sigma^2)= N(0,1)
	//static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。

	call_times++;
	if (call_times % 2 == 1) {
		real_t r = 10;

		{
			x_i = (Multiplier * x_i) % Meisen;
			random1 = ((real_t)x_i / Meisen);  //现在d是0-1的分布。最小标准随机数生成器
			x_i = (Multiplier * x_i) % Meisen;
			random2 = ((real_t)x_i / Meisen);  //现在d是0-1的分布
			random2 *= PI_mul_2;
			r = sqrtf(-2 * logf(random1)); //刷新r
		}

		gauss_1 = r * cosf(random2);
		gauss_2 = r * sinf(random2);    //一次可以计算两个高斯随机变量

		return gauss_1 * sigma + ave;
	}
	else
	{
		return gauss_2 * sigma + ave;
	}
}



//#define Meisen 2147483647      //梅森素数（形式为2^(p)-1的素数）：m=2^(31)-1
//#define Multiplier 16807       //乘子 a=7^(5)

/// <summary>
/// 高斯随机变量生成函数。使用Box_Muller方法，避免三角函数运算，速度稍微快一些
/// </summary>
/// <returns>高斯随机数，均值0，方差1</returns>
real_t random_gauss()
{
	//声明静态局部变量，这个变量实际地址不在函数中，程序编译时旧已经赋值了，调用这个程序的时候不会重复赋值。
	static unsigned int call_times = 0;//随机数调用次数。
	static real_t random1 = 0;
	static real_t random2 = 0; //0-1分布的随机数
	static real_t gauss_1 = 0;
	static real_t gauss_2 = 0; //均值0，方差1的高斯分布。 N(u,sigma^2)= N(0,1)
	static unsigned long int x_i = 271828;  //0-m之间的整数。是余数。  这个初始值是种子。

	call_times++;
	if (call_times % 2 == 1) {
		real_t r = 10;
		while (r > 1) // Box_Muller方法是一种拒绝方法。r如果大于1就淘汰
		{
			x_i = (Multiplier * x_i) % Meisen;
			random1 = ((real_t)x_i / Meisen);  //现在d是0-1的分布。最小标准随机数生成器
			x_i = (Multiplier * x_i) % Meisen;
			random2 = ((real_t)x_i / Meisen);  //现在d是0-1的分布

			random1 = random1 * 2 - 1;   //将0 - 1分布 变成-1 - 1的分布
			random2 = random2 * 2 - 1;


			r = random1 * random1 + random2 * random2; //刷新r
		}
		r = sqrtf(-2 * logf(r) / r);
		gauss_1 = r * random1;
		gauss_2 = r * random2;    //一次可以计算两个高斯随机变量

		return gauss_1;
	}
	else
	{
		return gauss_2;
	}
}



